Method of earthquake prediction

ABSTRACT

A method of predicting an earthquake in a seismically active region below an inosphere, including the steps of measuring ( 10, 12 ) a relative fluctuation of plasma density ( 26 ) in the inosphere ( 20 ), inferring a relative amplitude of an acousto-gravity wave in the inosphere ( 20 ), and inferring an earthquake magnitude from the relative amplitude of the acousto-gravity wave.

FIELD AND BACKGROUND OF THE INVENTION

The present invention relates to a method for predicting earthquakes.

Large earthquakes are among the most destructive of natural disasters.The toll that these disasters have taken in both lives and property iswell known and need not be recited here. Many methods have been proposedin the past for predicting earthquakes, to provide enough warning forthe people affected to prepare. One geophysical method, for example,based on the gravitational field turbulence associated with a seismicevent, was proposed by Chiba (Jiro Chiba, Proceedings of The Instituteof Electrical And Electronic Engineers, 1993 International CarnahanConference on Security Technology, pp. 215-218). Another geophysicalmethod, based on measurements of natural low frequency radio signals,was proposed by Hayakawa (M. Hayakawa, Phys. Earth. and Plan. Int., vol.77 pp. 127-135 (1993)). These methods provide only about 20 secondsadvance notice of an earthquake.

There is thus a widely recognized need for, and it would be highlyadvantageous to have, a method for predicting earthquakes that providesmore advance warning than known methods.

SUMMARY OF THE INVENTION

According to the present invention there is provided a method ofpredicting an earthquake in a seismically active region below anionosphere, comprising the steps of: (a) measuring a relativefluctuation of plasma density in the ionosphere; (b) inferring arelative amplitude of an acoustic-gravity wave in the ionosphere; and(c) inferring an earthquake magnitude from the relative amplitude of theacoustic-gravity wave.

The present invention is based on the discovery of an earthquakeprecursor that appears in the ionosphere above the epicenter severalhours before the earthquake. FIGS. 1A through 1F show two way traveltimes of reflections from the ionosphere, as a function of frequency,measured by a vertical ionosonde, in the 24 hours before an earthquakethat occurred in the Vrancha region of Romania on May 30-31, 1990. FIG.1A shows the structure of the normal ionosphere, 24 hours before theearthquake. Only one reflecting layer is present, and it reflects onlyat frequencies up to about 15 MHz. As the time of the earthquakeapproaches, more reflecting layers develop, and they reflect over awider range of frequencies. Because the reflection frequency isproportional to the square root of the ion density, this wider range offrequencies indicates that the ionosphere becomes more inhomogeneousabove the epicenter in the 12 to 16 hours before the earthquake. Theseinhomogeneities are believed to be produced by acoustic-gravity waves(AGW) associated with the buildup of strain in the Earth's crustimmediately before the earthquake.

This phenomenon provides up to 20 hours warning in advance of anearthquake, but it does not provide a measure of the magnitude of theimpending earthquake. According to the present invention, a measure ofthe magnitude of the earthquake is provided by measurements of obliquescattering of radio waves from the ionospheric inhomogeneities. Thesemeasurements are interpreted using a theory of the coupling of AGW andionospheric plasma density that is presented herein.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention is herein described, by way of example only, withreference to the accompanying drawings, wherein:

FIGS. 1A through 1F show reflections from the ionosphere, as a functionof frequency, at 2 hour decrements in the 24 hours before an earthquake;

FIG. 2 is a schematic illustration of the measurement geometry of thepresent invention;

FIG. 3 is a flow chart of the present invention.

MODULATION OF PLASMA DENSITY BY ACOUSTIC-GRAVITY WAVES

During its movement through the ionospheric cold dense plasma, an AGW offrequency ω and wavenumber vector {right arrow over (k)} causes theneutral particle density N_(m) and the rate of ionization per unitvolume q({right arrow over (r)}) to vary as follows:

N _(m) =N _(m0) +δN _(m) exp[i({right arrow over (k)}·{right arrow over(r)}−ωt)]  (1)

q({right arrow over (r)})=q ₀ +δq exp[i({right arrow over (k)}·{rightarrow over (r)}−ωt)]  (2)

(A. Gossard and Y. Huk, Waves in the Atmosphere, 1975). Here, δN_(m) andδq are the disturbances of neutral particle density and of the rate ofionization per unit volume due to the AGW, {right arrow over (r)} is theposition vector, and t is time. Let us now estimate the changes inplasma concentration caused by the AGW. Here we must consider two cases,the electrostatic case and the electrodynamic case.

In the D-layer and lower E-layer, the problem is electrostatic. At thesealtitudes, the ionization rate is as in equation (2), but therecombination rate is of the form α_(ef)N_(e) ², where α_(ef) is therecombination coefficient, on the order of 0.2×10⁻⁶ cm³s⁻¹ to 10⁻⁶cm³s⁻¹. At these altitudes (50-100 km) ions move with the speed ofneutral particles, V^(i)≈V^(m)≈50-100 m/s, and the speed of chargedparticles is less than the speed of the AGW. In other words, the chargedparticle motion has no influence on the change of plasma concentration,i.e., V^(i)<<V_(A), where V_(A) is the speed of sound. In this case,from the equation of conservation of particles, $\begin{matrix}{\frac{\partial N}{\partial t} = {q - {\alpha_{ef}N^{2}} - {{div}( {N{\overset{arrow}{V}}_{\alpha}} )}}} & (3)\end{matrix}$

it follows that the characteristic time of change of particle densitydue to recombination, τ_(R), is much less than the characteristic timeof particle transport, τ_(T). In fact, for N_(e0)=N_(i0)=N₀=10¹⁰m⁻³,V_(m)≈50 m/s and the characteristic scale of quasiregular plasmainhomogeneity Λ≈200 km at these altitudes, we obtainτ_(R)=(α_(ef)N)⁻¹≈100 s-500 s and τ_(T)=Λ/V^(m)≈4000 s. Withτ_(T)<<τ_(R), we can ignore the movement of charged particles in (3),i.e., $\begin{matrix}{\frac{\partial N}{\partial t} = {q - {\alpha_{ef}N^{2}}}} & \text{(3a)}\end{matrix}$

Introducing (2) to (3a) and linearizing, we find that the amplitude ofplasma inhomogeneities due to AGW can be estimated as: $\begin{matrix}{{\delta \quad N_{e,i}} = {{N_{e,i} - N_{{0e},i}} \approx {\delta \quad q2\alpha_{ef}\frac{N_{{0e},i}}{\omega^{2} + ( {2\alpha_{ef}N_{{0e},i}} )^{2}}}}} & (4)\end{matrix}$

Here, q₀=α_(ef)N_(0e) ² (N_(0e)≈N_(0i), plasma is quasi-neutral),δq=q₀δN_(m)/N_(m0)=q₀δ_(A), where δ_(A) is the relative amplitude of theAGW. Using this notation, we finally have: $\begin{matrix}{{\delta \quad N_{e,i}} = {\delta_{A}2\alpha_{ef}\frac{N_{{0e},i}^{3}}{\omega^{2} + ( {2\alpha_{ef}N_{{0e},i}} )^{2}}}} & (5)\end{matrix}$

for AGW with frequency ω<2α_(ef)N₀, δN/N₀=δ_(A)/2. Hence, for lowerfrequencies, the plasma disturbances are proportional to their source,i.e., to the amplitude of the AGW. For higher frequencies, i.e.,ω>2αefN₀, δN/N₀=2δ_(A)F(ω), where F(ω)=(α_(ef)N₀/ω)². Because in theD-layer and the lower E-layer, α_(ef)N₀=0.002 s⁻¹-0.01 s<<ω, F(ω)<<1,and δN/N₀<<δ_(A). Hence, for high frequencies, plasma disturbances areweak, i.e., much less than the amplitude of the AGW.

In the F-layer, the problem is electrodynamic. For these altitudes, theamplitude of the moving plasma inhomogeneities depends on theorientation of the wave vector {right arrow over (k)} of the AGWrelative to the geomagnetic field {right arrow over (B)}₀ and relativeto the plasma drift velocity {right arrow over (V)}_(d). As was shown inC. O. Hines, Can. J. Phys. Vol. 38 pp. 1441-1481 (1960), the amplitudeof moving plasma disturbances (MPD) is maximal when the condition ofspatial resonance is “working” (i.e., the plasma drift velocity is equalto the AGW phase velocity):

ω={right arrow over (k)}·{right arrow over (V)} _(d)  (6)

The MPD are created by the transport processes during interactions ofcharged particles with neutrals. The latter are modulated according toequation (1) by the AGW, so the charged particles also are modulated bythe AGW. Because at these altitudes (>150 km) magnetic and electricfields are significant, the interaction between the AGW and the plasmais electrodynamic. In this case, because the frequency of ion-neutralinteractions, ν_(im), is on the order of N_(m)(T_(m)+T_(i))^(½), whereT_(m) is the temperature of the neutral species and T_(i) is thetemperature of the ions, the modulation of N_(m) and T_(m) by AGWproduces a corresponding change in ν_(im) and hence the correspondingchange in the ion velocity {right arrow over (V)}_(i). This modulationof {right arrow over (V)}_(i) causes the redistribution of plasmadensity and the creation of δN.

Let us consider that the plasma is isothermal, i.e., T_(e)=T_(i)=T,which is the case in the lower ionosphere, up to about 200 km, andchoose a coordinate system in which the magnetic field {right arrow over(B)}₀ is parallel to the z-axis. The particle movement is considered ina coordinate system (x,y,z) which is at rest with respect to averageneutral flow. The background plasma is quasi-regular and homogeneous. If{right arrow over (V)}₁ ^(m) is the velocity of neutral species in thefield of the AGW, {right arrow over (k)} is the wavevector of the AGW,δ_(A)=δN_(m)/N_(m0)=V₁ ^(m)/V_(ph), where V_(ph)=ω/k is the phasevelocity of the AGW, and δ_(T)=δT_(m)/T_(m0) is the relative amplitudeof the modulation of the temperature by the AGW, then we can accuratelyrepresent the frequency of interaction as a periodic function due tomodulation cause by the AGW (D. F. Martin, Proc. Roy. Soc. Canada, Vol.A209, p. 216 (1950)):

ν_(im)=ν_(im) ⁰[1+(δ_(A)+δ_(T)/2)exp{i({right arrow over (k)}·{rightarrow over (r)}−ωt)}]=ν_(im) ⁰+δν_(im)  (7a)

ν_(ei)=ν_(ei) ⁰[1+(δ_(i)+3δ_(T)/2)exp{i({right arrow over (k)}·{rightarrow over (r)}−ωt)}]=ν_(ei) ⁰+δν_(ei)  (7a)

Because of the difference between {right arrow over (V)}^(e) and {rightarrow over (V)}^(i) in the magnetic ({right arrow over (B)}₀) andelectric ({right arrow over (E)}₀) fields (E₀≈10 mV/m-50 mV/m) betweenparticles (electrons and ions), the ambipolar field is created, i.e.,δ{right arrow over (E)}=−∇φ (φ is the potential of the ambipolar field),and also the changes of the tensors of particle mobility:

{circumflex over (M)} ^(α) ={circumflex over (M)} ₀ ^(α) +δ{circumflexover (M)} ^(α) ,α=e,i

{right arrow over (E)}={right arrow over (E)} ₀ +δ{right arrow over(E)}  (8)

where $\begin{matrix}{{\delta \quad {\overset{arrow}{M}}^{\alpha}} = {\frac{\delta \quad \nu_{\alpha}}{\nu_{\alpha}^{0}}\begin{pmatrix}M_{p}^{\alpha} & 0 & 0 \\0 & M_{p}^{\alpha} & 0 \\0 & 0 & {- M_{\parallel}^{\alpha}}\end{pmatrix}}} & (9)\end{matrix}$

Here, p denotes the Pedersen component of electrons (α=e) and ions (α=i)mobility, i.e., in the direction of the electric field, and M_(∥) ^(α)is the mobility of charged particles along the magnetic field. Finally,in the electrodynamic case, we can present (taking into account that alldisturbances are weaker than the background values) in the linear case,the disturbances of charged particle (electrons and ions) velocities (M.G. Gel'berg, Geom. And Aeron. Vol. 20, pp. 271-274 (1980)).$\begin{matrix}{{\delta \quad {\overset{arrow}{V}}^{\alpha}} = {{q( {{\delta \quad {\hat{M}}^{\alpha}{\overset{arrow}{E}}_{0}} + {{\hat{M}}_{0}^{\alpha}\delta \quad \overset{arrow}{E}}} )} + {\frac{B_{0}}{Q_{\alpha}}{\hat{M}}_{0}^{\alpha}\delta \quad {\overset{arrow}{V}}^{m}} - {\frac{D_{\alpha}B_{0}}{Q_{\alpha}}{\hat{M}}_{0}^{\alpha}{\nabla\delta}\quad \frac{N}{N_{0}}}}} & (10)\end{matrix}$

Here, D_(α) is the coefficient of diffusion of electrons (α=e) or ions(α=i); Q_(α)=ω_(Hα)/ν_(α), ω_(Hα) is the gyrofrequency if plasmaparticles in the magnetic field; ν_(e)=ν_(em)+ν_(ei);ν_(i)=ν_(im)+ν_(ei). In the case when disturbances of the plasmaparticles due to their modulation by the AGW are weak, i.e., δN<<N₀, wecan in equation (3) exclude the ambipolar field. Linearizing equation(3) and taking into account equation (10), an equation can be obtainedwhich, for longitudinal AGW, gives $\begin{matrix}{\frac{\delta \quad N}{N_{0}} = {( {\omega - {\overset{\_}{k} \cdot {\overset{\_}{V}}_{d}}} )\frac{{\delta_{A}{kV}_{ph}\cos^{2}\chi} + {( {\delta_{A} + {\delta_{T}/2}} )( \frac{\overset{arrow}{k} \cdot {\overset{arrow}{E}}_{0}}{Q_{i}B_{o}} )}}{( {\omega - {\overset{arrow}{k} \cdot {\overset{arrow}{V}}_{d}}} )^{2} + ( {D_{\alpha}k_{\parallel}^{2}} )}}} & (11)\end{matrix}$

Here, {right arrow over (V)}_(d) is the drift velocity of chargedparticles, and χ is the angle between {right arrow over (k)} and {rightarrow over (B)}₀. From (11) it follows that ifQ_(i)V_(ph)/V_(d)cos²χ<<1, then the influence of modulation exceeds theeffects of interaction between the plasma and the neutral species. Inthis case, the angle ψ={fraction (π/2+L )}−χ is less thanψ_(k)=sin⁻¹{square root over (V_(d)+L /Q_(i)+L V_(ph)+L )}. For analtitude of 200 km, Q_(i)=10⁻¹ to 10⁻², V_(ph)˜600 ms⁻¹, E₀=30 mV/m, andwe obtain ψ_(k)˜7° to 20°. Expression (11) has a maximum for AGWsatisfying the condition ω−{right arrow over (k)}·{right arrow over(V)}_(d)≈D_(α)k_(∥) ², where k_(∥) is the component of {right arrow over(k)} parallel to the magnetic field: $\begin{matrix}{( \frac{\delta \quad N}{N_{0}} )_{\max} = \frac{{\delta_{A}k\quad V_{p\quad h}\cos^{2}\chi} + {( {\delta_{A} + {\delta_{T}/2}} )( {{\overset{arrow}{k} \cdot {{\overset{arrow}{E}}_{0}/Q_{i}}}B_{0}} )}}{D_{\alpha}k_{\parallel}^{2}}} & (12)\end{matrix}$

The above analyses illustrate the principal possibility ofelectrodynamic and electrostatic mechanisms of the redistribution ofplasma by longitudinal AGW and, finally, of the creation of plasmainhomogeneities.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention is of a method of predicting earthquakes as muchas 12-16 hours in advance.

The principles and operation of earthquake prediction according to thepresent invention may be better understood with reference to thedrawings and the accompanying description.

The present invention is based on N. D. Philipp and N. Sh. Blaunshtein,Power of H_(E)-Scatter Signals, Proc. XII Union Conference on RadiowavePropagation, Kharkov, USSR, 1978, pp. 166-169, which is incorporated byreference for all purposes as if fully set forth herein.

Referring now to the drawings, FIG. 2 illustrates schematically thegeometric setup of the apparati used in the present invention. A radiofrequency transmitter 10 and a radio frequency receiver 12 are separatedby up to 1000 km or more on the surface of the Earth. Together,transmitter 10 and receiver 12 constitute an oblique ionosonde.Transmitter 10 transmits radio frequency energy in a solid angle 22.Receiver 12 receives radio frequency energy from a solid angle 24. Inparticular, receiver 12 receives radio frequency energy from transmitter10 that is reflected from anomalous plasma inhomogeneities in theportion of an ionospheric anomaly 20 that lies within the volume ofoverlap 26 of solid angles 22 and 24. The separation of transmitter 10and receiver 12 is chosen to be consistent with the elevations of solidangles 22 and 24, the power of transmitter 10 and the sensitivity ofreceiver 12; the preferred separations are between 500 km and 1000 km.The radio frequency energy should be transmitted at a frequency greaterthan the ionospheric plasma frequencies, so frequencies greater thanabout 20 MHz are used.

Between transmitter 10 and receiver 12 are deployed one or more verticalionosondes 14. When an anomalous reflection pattern such as thatillustrated in FIG. 1D is detected by one of vertical ionosondes 14,transmitter 10 and receiver 12 are activated to measure the relativefluctuation of plasma density, δN/N₀, of equation (11), using equation(6) of Philipp and Blaunshtein: $\begin{matrix}{P_{R} = {P_{T}\frac{{{\alpha\lambda}^{3}( {\delta \quad {N/N_{0}}} )}^{2}}{32e\sqrt{2\pi}r^{4}}{\int_{V}{\frac{G_{T}G_{R}\sin^{2}}{\sin^{2}( {\theta/2} )}{\exp ( {{- 8}\pi^{2}\alpha^{2}\psi^{2}\sin^{2}\frac{\theta}{2}} )}{V}}}}} & (13)\end{matrix}$

In this equation, P_(T) is the transmitter power, P_(R) is the receiverpower, r is the distance from transmitter 10 and receiver 12 to overlapvolume 26 (if transmitter 10 and receiver 20 are not equidistant fromoverlap volume 26, the term r⁴ should be replaced with r_(T) ²r_(R) ²,where r_(T) is the distance from transmitter 10 to overlap volume 26 andr_(R) is the distance from receiver 12 to overlap volume 26), λ is thewavelength of the radio frequency radiation, G_(T) is the transmittergain, G_(R) is the receiver gain, the angles θ, ψ and χ are defined inPhilipp and Blaunshtein, and the integration is performed numericallyover overlap volume 26. (The relative fluctuation of plasma density isdenoted as δN/N₀ herein, instead of the notation “ΔN/N” used in Philippand Blaunshtein.) Note that the latest arrival time measured by verticalionosonde 14 provides an upper bound on the altitude of anomaly 20, andthe integration in equation (13) should be limited to altitudes belowthis upper bound. α is the ratio L/λ, where L is the longitudinal extent(along the earth's magnetic field) of the inhomogeneities in anomaly 20.This ratio may be estimated using equation (11) of Philipp andBlaunshtein, as described therein.

Solid angles 22 and 24 are scanned to provide a map of δN/N₀ in theionosphere between transmitter 10 and receiver 12. The epicenter of theimpending earthquake is predicted to be directly below the point ofmaximum δN/N₀, to within about 5 km. δ_(A) is estimated using equation(12), on the assumption that δ_(T) is much smaller than δ_(A). Theionospheric parameters of equation (12), such as wave vector {rightarrow over (k)}, are well known to those ordinarily skilled in the art,and may be found, for example, in M. G. Gel'berg, The Inhomogeneities ofthe High-Latitude Ionosphere, Nauka, Novosibirsk, 1986. To ensure thevalidity of an analysis based on equation (12), the integration inequation (13) is limited to altitudes at which the electrodynamic modelis valid, i.e., above about 80 km. It has been found empirically, forlow-magnitude earthquakes, that for earthquakes of magnitude greaterthan about 3 on the Richter scale, the earthquake magnitude isapproximately 0.8 times δ_(A). If the predicted magnitude is greaterthan about 5, the populace near the predicted epicenter is warned toevacuate or take other precautionary measures.

Earthquakes tend to occur in fault zones, which are linear geologicalfeatures. Therefore, transmitter 10 and receiver 12 preferably aredeployed at opposite ends of a fault zone, and are used to probe theionosphere above the fault zone. For example, to monitor the San Andreasfault of California, transmitter 10 may be deployed in San Diego andreceiver 12 may be deployed in San Francisco (or vice versa). In thecase of a relatively restricted area of seismic activity, such as theVrancha Region of Romania, transmitter 10 and receiver 12 may bedeployed on opposite sides of the area of seismic activity. In thelatter case, it is not necessary to scan solid angles 22 and 24.Instead, solid angles 22 and 24 are kept fixed, and are made wide enoughfor overlap volume 26 to cover the entire area of seismic activity.

As an alternative to deploying an oblique ionosonde, a sufficientlydense array of vertical ionosondes 14 may be deployed, with eachvertical ionosonde 14, in addition to monitoring the reflectivitystructure of the ionosphere thereabove, also serving as a combinedtransmitter and receiver to measure δN/N₀ in a solid angle directlythereabove. The upper limit of integration for equation (13) is providedby the highest altitude inferred from the maximum two way travel time ofthe reflection patterns received by vertical ionosondes 14. The lowerlimit of integration for equation (13) is provided by either the lowestaltitude inferred from the minimum two way travel time of the reflectionpatterns received by vertical ionosondes 14, or by the lowest altitudeof validity of the electrodynamic model (˜80 km), which ever is higher.

FIG. 3 is a flow chart of the present invention. In box 100, verticalionosondes 14 monitor the reflectivity structure of the ionosphere foranomalies. Upon detection of an anomalous reflection pattern (box 102),transmitter 10 and Receiver 12 are activated to measure δN/N₀ (box 104).δ_(A) is estimated from δN/N₀ in box 106, and the magnitude of theearthquake is estimated in box 108.

While the invention has been described with respect to a limited numberof embodiments, it will be appreciated that many variations,modifications and other applications of the invention may be made.

What is claimed is:
 1. A method of predicting an earthquake in aseismically active region below an ionosphere, comprising the steps of:(a) measuring a relative fluctuation of plasma density in the ionosphereprior to the earthquake; (b) estimating a relative amplitude of anacoustic-gravity wave in the ionosphere from said measured relativefluctuation; and (c) estimating an earthquake magnitude from saidrelative amplitude of said acoustic-gravity wave.
 2. The method of claim1, wherein said measuring of said relative fluctuation of plasma densityis effected by scattering radio frequency energy from the ionosphere. 3.The method of claim 2, wherein said radio frequency energy has afrequency greater than about 20 MHz.
 4. The method of claim 2, whereinsaid measuring is effected using a transmitter that transmits said radiofrequency energy obliquely in a first solid angle towards the ionosphereand a receiver that receives radio frequency energy obliquely in asecond solid angle from the ionosphere, said relative fluctuation ofplasma density being measured in an overlap volume defined by said firstsolid angle and said second solid angle.
 5. The method of claim 4,wherein said transmitter and said receiver are separated by betweenabout 500 km and about 1000 km.
 6. The method of claim 4, wherein saidmeasuring includes a scanning of at least one solid angle selected fromthe set consisting of said first solid angle and said second solidangle, the method further comprising the step of: (d) creating a map ofsaid relative fluctuation of plasma density, said map having a point ofmaximum relative fluctuation of plasma density, the earthquake beingpredicted to occur substantially directly below said point of maximumrelative fluctuation of plasma density.
 7. The method of claim 4,wherein said transmitter and said receiver are deployed at opposite endsof the seismically active region.
 8. The method of claim 4, wherein saidtransmitter and said receiver are deployed outside the seismicallyactive region, with the seismically active region lying between saidtransmitter and said receiver.
 9. The method of claim 1, wherein saidmeasuring of said relative fluctuation of plasma density is effected ina plurality of volumes of the ionosphere, the method further comprisingthe step of: (d) creating a map of said relative fluctuation of plasmadensity, said map having a point of maximum relative fluctuation ofplasma density, the earthquake being predicted to occur substantiallydirectly below said point of maximum relative fluctuation of plasmadensity.
 10. The method of claim 9, wherein said measuring of saidrelative fluctuation of plasma density is effected by scattering radiofrequency energy from said volumes, said radio frequency energy beingemitted by a transmitter in a first solid angle and being received by areceiver in a second solid angle, said plurality of volumes beingdefined by said first solid angle and said second solid angle as atleast one solid angle selected from the set consisting of said firstsolid angle and said second solid angle is scanned.
 11. The method ofclaim 10, wherein said radio frequency energy has a frequency greaterthan about 20 MHz.
 12. The method of claim 10, wherein said transmitterand said receiver are separated by between about 500 km and about 1000km.
 13. The method of claim 1, further comprising the step of: (d)monitoring the ionosphere for an anomalous reflection pattern; saidmeasuring of said relative fluctuation of plasma density being initiatedupon detection of said anomalous reflection pattern.
 14. The method ofclaim 13, wherein said monitoring is effected using a verticalionosonde.
 15. The method of claim 14, wherein said measuring of saidrelative fluctuation of plasma density is effected by scattering radiofrequency energy from said ionosphere, using a transmitter thattransmits said radio frequency energy obliquely in a first solid angletowards the ionosphere and a receiver that receives radio frequencyenergy obliquely in a second solid angle from the ionosphere, saidvertical ionosonde being deployed substantially between said transmitterand said receiver.
 16. The method of claim 15, wherein said radiofrequency energy has a frequency greater than about 20 MHz.
 17. Themethod of claim 15, wherein said transmitter and said receiver areseparated by between about 500 km and about 1000 km.
 18. The method ofclaim 13, wherein said monitoring is effected using a plurality ofvertical ionosondes.
 19. The method of claim 18, wherein each of saidreflection patterns detected by each of said plurality of ionosondes hasa minimum two way travel time defining a minimum altitude and a maximumtwo way travel time defining a maximum altitude, and wherein saidmeasuring of said relative fluctuation of plasma density is effected bytransmitting radio frequency energy vertically in a solid angle fromeach of said plurality of ionosondes and receiving said radio frequencyin said solid angle by said each of said plurality of ionosondes; foreach of said plurality of ionosondes, said solid angle, said maximumaltitude, and the greater of said minimum altitude and a lowest altitudeof validity of an electrodynamic model defining a volume in theionosphere, the method further comprising the step of: (e) creating,from said measurements of said relative fluctuation of plasma density insaid volumes, a map of said relative fluctuation of plasma density, saidmap having a point of maximum relative fluctuation of plasma density,the earthquake being predicted to occur substantially directly belowsaid point of maximum relative fluctuation of plasma density.
 20. Themethod of claim 1, wherein said measuring is effected at least about 12hours before the earthquake.